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ABSTRACT 


Recent developments In the numerical simulation of turbulent flows are discussed. 
This limited survey covers computational requirements for the direct simulation of tur- 
bulence » simulation of arbitrary homogeneous flows, a new expansion technique for wall- 
bounded flows with application to pipe flow, and possibilities of flow representations 
or modeling techniques that allow the simulation of high-Reynolds-number flows with a 
relatively small number of dependent variables. 


INTRODUCTION 


The numerical simulation of turbulent flows has a short history. Around 1950 
von Neumann [1] and Emmons [2] proposed an attack on the turbulence problem by numerical 
simulation. But one could point to a beginning 20 years later in 1970-71 when Deardorff 
[3] reported on a large-eddy simulation of -turbulent channel flow on a 24 x 20 x 14 mesh 
and a direct simulation of homogeneous, isotropic turbulence was accomplished on a 
32^ mesh by Orszag and Patterson [4], Perhaps the arrival of the CDC 6600 triggered 
these initial efforts. Since that time a number of developments have occurred along sev- 
eral fronts. Of course, faster computers with more memory continue to become available 
and it appears that this trend will not let up at least for another 3-5 years. In addi- 
tion, new algorithms have been developed which extend or improve capabilities in turbu- 
lence simulation. For example, the simulation of arbitrary homogeneous flows and the 
efficient simulation of wall-bounded flows are now possible. Finally, the prospects for 
simulating high-Reynolds-nimiber flows with limited computational resources have been 
realized with the development of sub grid -mode ling techniques and vortex methods. As a 
result of all these developments, applications of turbulence simulation have been increas- 
ing rapidly in the past few years and will continue to do so. 

In this paper I present examples of these new developments and discuss prospects for 
future developments. 


PROBLEM OF NUMERICAL SIMULATION 


We consider an incompressible flow whose time evolution is given by the Navier-Stokes 
equations for the velocity, u(x,t), and the pressure, p(x,t), as 

3u 2 

+ u • 7u « -7p + v7^u , (la) 

at—— — 

V • u « 0 , (lb) 

along with appropriate initial and boundary conditions. It is assumed that the 
density ■ 1. The character of the solution depends on the Reynolds number of the flow, 

Re « UL/v, where U and L are a characteristic velocity and length of the large scale 
and V is the kinematic viscosity. For small Reynolds numbers, one obtains a laminar 
flow that is smoothly varying in space and time; for large Reynolds numbers, one obtains 
a turbulent flow. Turbulent flows have been described as random, chaotic, vortical, 
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three-dimensional, and unsteady, and they are known to contain a wide range of scales* 

It is the combination of all these attributes that makes the numerical simulation of such 
flows extremely challenging. 

In turbulent pipe flow, for example, we estimate, according to universal equilibrium 
theory [5], the smallest important scale of turbulence to be proportional to the dissipa- 
tion, or Kolmogorov, length, n • where e is the energy-dissipation rate per 

unit mass, and the largest important scale to be some multiple of the pipe diameter. 

Using the volume-averaged e given by 


e - - 


4Uv au 
D 3r 


wall 


where U is the mean velocity and D is the pipe diameter, we find that 

n - 

where Re - UD/v and f is the friction factor, 

f - 8 (u,/U)2 , 


and is the wall shear velocity given by 


2 3U 


wall 


The friction factor, given implicitly by the formula [6], 



is only weakly dependent on Reynolds number so that the required number of mesh points 
on a three-dimensional grid would be proportional to (D/n)^ « Re^/**. Figure 1 shows the 
energy spectrum measurements of Laufer [7] for high-Reynolds-number (Re * 500,000) pipe 
flow. The pipe diameter is 25.4 cm. The wave number corresponding to the Kolmogorov 
length, k^ ■ 2Tr/n, is seen to be well beyond the measured data. To simulate reliably the 
dissipation of turbulence energy, the grid spacing must be somewhat smaller than the \ 
length scale corresponding to the peak in the dissipation spectrum. If isotropy of the 
small scales is assumed the dissipation spectrum is proportional to k^EiCkj). In 
Laufer' s experiment this peak, away from the wall, corresponds to a length of 150n or 
O.bSD, Lawn’s experiments [8] with lower Reynolds numbers (Re ■ 37,000 to 250,000) and 
those of Bakewell and Lumley [9] (Re ■ 8,700) indicate that the location of the peak 
dissipation scales on diameter (l.e., occurs at ->0.03D) , not on Kolmogorov length. (It 
would appear that these Reynolds numbers are coo low for universal equilibrium theory to 
apply. Nevertheless, I will assume that the Kolmogorov length is still a good measure 
of the smallest turbulent scale.) Thus, at low Reynolds numbers, the Kolmogorov wave 
number is much closer to the peak of the dissipation spectrum. For example, if 
Re - 5,000 the dissipation peak would correspond to a length of »7n. 

Therefore, as an estimate of the mean spacing between grid points. A, required in 
the direct simulation of turbulent pipe flow, we take A « 2n. Table 1 gives correspond- 
ing estimates of the number of mesh points required for several Reynolds numbers, assum- 
ing that Che computational domain extends 10 diameters in the streamwlse direction. 

(This estimate could be off by a factor of 3 either way. Some measurements and their 
interpretation suggest correlation lengths of 20D, others correlation lengths of 2D; 
see Coles [10].) It appears that only the lowest Reynolds number case is accessible ' to 
present day supercomputers (10® floating point operations per second, 10^ words storage). 

In addition it should be verified that the spacing A « 2n is sufficiently small to 
allow resolution of all important turbulence phenomena near the wall. The grid spacing 
measured in wall units is given by 
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For Re from 5,000 Co 500,000, ranges from 3.2 to 7.6 (see table 1). This spacing 
should provide sufficient resolution to reproduce all important wall- layer structures 
(such as streamwise streaks), which have characteristic lengths of 50-100 wall units. 

See Moin [H] further discussion of this point. 

The number of time-steps, Ns, required to follow one realization for a time T and 
obtain reasonable statistics also depends on Reynolds number. The time-step At is 
roughly limited to 


At 



Using the above estimate for A 


and (2/f)^/** 5 3 


we find that 


At 


< ---- . 
URe^^“ 


And if T - 100 D/U, then 


Ns - ^ > URe^/" , 


or 10,000 steps for Re • 5,000. 

DIRECT SIMULATION OF TURBULENCE 


Homogeneous flows ; 

Probably the most ambitious direct simulations of turbulent flows to dace are the 
computations of a variety of homogeneous turbulent flows by Rogallo [12] on a 128^ mesh. 
Rogallo writes the velocity field u as the sum of a mean component and a turbulent 
component , 


u ■ U + u* , 

and assumes that the components of U have the form 

“m ■ “mn(t)Xn . 

where repeated indices are summed. By transforming to coordinates, x, moving with the 
mean flow he obtains momentum and continuity equations for u* that contain no explicit 
dependence on x. Then Rogallo assumes chat the turbulent component of the homogeneous 
flow is periodic in ic-space, with period in direction m; therefore, no further 

boundary conditions are required, and spatial derivatives can be computed accurately by 
Fourier interpolation. Thus, u*(x,t) is represented by 


u*(x,t) « u(k,t)e^-“- . 

k ~ “ 

Here the mth component of k is 

and I ranges over -N/2 + L<i<N/2-l, where N « 128 for Rogallo *s simulations 
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and the Navler-Stokes equations become a 3(N - 1)^ system of ordinary differential 
equations (ODEs)t 







- vk'Sni • 


where _ k* - kjjkjj , and 


5V^(k) - Z Uj,(k')Ujj(k - k') . 
k* 

(In the above, zero mean flow, Uffln - 0. has been assumed in order to simplify the pres- 
entation.) To avoid explicit evaluation of the convolution sums 1%% requiring 0(N®) 
operations per step, fast Fourier transforms (FFTs) are used to return to physical space 
where the required products are formed and then transformed (by FFTs) back to Fourier 
space. Consequently, only 0(N^ log N) operations per step are reqiiired. 


Suppose the tensor ^nn is decomposed into a symmetric (strain, R) and antisymmet- 
ric (vorticity, f^) tensors; then the only constraint on R is that it be traceless 
R^j^ ■ 0, but the vorticity must satisfy the evolution equation 

\n ^nk^kn ^k^kn “ ° • 


Rogallo simulated four flows: plain strain, U22 “ *^33 “ const; axisymmetric strain, 

U22 * U33 ■ -(1/2)1)22 “ const; shear, U22 • const; and rotation, U23 • “U31 ■ const; 

Figure 2 shows the energy spectrum and the component spectra for the case of homo- 
geneous shear, s U22 several dimensionless times, st, after starting from an iso- 
tropic field at St " 0. Note that at the .small scales (large k) the component energies 
are nearly equally distributed, whereas the*streamwise component (E^) highest 

large-scale energy, and the shear-direction component (Ej) has the least. Figure 3 shows 
the tendency of the nondimensional turbulence energy to decrease and then equilibrate. 

All these features of the simulation are in agreement with experiment. 

Agreement with experiment is of course desirable and provides confidence in the 
solutions but is not the ultimate goal of these simulations. New Information about tur- 
bulent flows or new insights into the nature of turbulence is the desired outcome.' For 
example, in the development of turbulence models, one hypothesizes relations between 
higher-order /multipoint moments of the velocity field and lower-order/single-point moments. 
Any desired statistics are readily computed from numerical simulations so that various 
hypotheses may be tested or new ones may suggest themselves. This is demonstrated in 
figure 4. The correlation between — the sum of the slow pressure strain term (cubic 
in turbulence velocities) and the deviator of the dissipation tensor and the energy 
tensor bij , given by 




’ 3 



is shown in figure 4a. The correlation, suggested by the work of Lumley and Newman [13], 
is seen to be very good for the simulations of homogeneous shear and represents a signif- 
icant improvement over earlier attempts to correlate the total pressure strain with the 
energy tensor. If is modeled as 


and 3 is adjusted separately for each field to minimize the square error of the fit, 
then one obtains very good agreement of the model with measured values of , as shown 
in figure 4b. In addition, Rogallo found that the adjusted values of B correlate well 
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with Reynolds number, Rl, as shown in figure 4c. (It can be shown analytically that 
3 -»• 2 as R^ -*• 0.) Thus, the model above for ^ with 3 - promising. 

Wall"bounded shear flows : 

The direct simulation of wall-bounded shear flows presents a new set of difficulties 
because of the presence of the wall. Part of the problem is due to the physical pro- 
cesses near the wall. The thin shear layer next to the wall is continually breaking up 
via three-dimensional (3-D) inertial instabilities resulting in a violent 3-D wrinkling 
of the vortex layer. In the spanwise direction, the scale of the breakup is of the order 
of the thickness of the layer or tens of wall units and, perhaps, somewhat more in the 
streamwise direction. 

Another part of the problem Is algorithmic and is due to the no-sllp boundary con- 
dition at the wall. One can no longer use Fourier series for the spectral expansion in 
the direction normal to the wall. Rather, to obtain rapid convergence independent of 
boundary constraints, one should employ global polynomials related to the eigenfunctions 
of a singular Sturm-Liouville problem [14]. Ghebyshev and Legendre polynomials are popu- 
lar choices but, as shown below in the case of pipe flow, other choices may be preferable 
because of special conditions of the problem at hand. Whatever the choice may be, this 
change in basis functions complicates the imposition of the no-slip condition and the 
satisfaction of the continuity constraint. These two conditions become global con- 
straints on the expansion which are generally difficult or costly to impose simulta- 
neously. By contrast, in the simulation of homogeneous flows using Fourier expansions 
in all three directions, the boundary conditions are built into the expansion and the 
divergence-free constraint is satisfied by a simple projection which is local in wave- 
number space. 

In the following, I shall describe a new technique for overcoming the algorithmic 
difficulties described above, at least for wall-bounded flows in simple geometries. The 
technique consists of expanding the velocity field in terms of a set of divergence-free 
vector functions satisfying the boundary conditions. First we write the Navier-Stokes 
equations in rotational form: 

^ ^ X u * -7^p + *1 u^^ + \>7^u , (2a) 

7 • u - 0 . (2b) 

Here j » 7 x u is the vorticity. The boundary condition at a wall is 

Other boundary conditions, such as periodicity or inflow-outflow, are imposed as 

appropriate. 

The role of the pressure in incompressible flows is to enforce the continuity con- 
dition. This may be expressed formally by recalling that an arbitrary vector field f 
may be uniquely decomposed into a sum of a divergence- free field satisfying tangency at 
the boundary and the gradient of a potential, 

f = jjj + 7(|) , 

7 • V » 0 • n|wall * ® * 

Let ^ be the projection operator that accomplishes this decomposition, that is, 

^f - ^ ; 

applying the projection operator ^ to equation (2a) we obtain 

3u , 

X u) + ^(v7^u) , 

eliminating the dynamic pressure [15]. The above equation is the starting point for the 
numerical scheme described below. 
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Vector expansion method ; 


We write u as the expansion 


N 

u(x,t) ■ 2 • (3) 

n«i 


where each satisfies 


^ - 0 , (4) 

and the homogeneous boundary condition on u. We need to derive a system of DDEs for the 
coefficients aj^(t)(n - 1,2, • . .,N). We do so by substituting the expansion (3) into 
(2a) and caking the inner product of the result with a set of weight vectors 
(m « 1,2, . . .,N) satisfying 


V 




0 


(5) 


and 


-m 


-’wall 


( 6 ) 


That is, 


J t (x) • [momentum eq.]dV * 0 , m - 1,2, • • . ,N • (7) 

V - 

If Che fora a conplece sec and N •» <•, chis operaclon is equlvalenc co applying Che 

projecclon operator because 

f d (x) . 7* dV - -f (7 . c )4, d7 + f ♦(?„ • n)dS 

Jv •'s ‘ 

- 0 

for an arbitrary scalar field An explicit representation of ^ may be given as 


■ S SAn<5„ • . 

where 

0 ^ • <C • C • 

^n ' -n -m 

The result of the computation given in (7) is the following system of ODEs; 

Aa + vBa ■ g , (8) 


where 



J^dV . 


(7 X 7 X J(j^)dV 


and gjQ is qtiadratic in the a^'s. 


(9) 

( 10 ) 
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( 11 ) 


8ni ■ -J Sm • (y “)dV . 

Thus, each evaluation of the requires the solution of a linear system and the 

computation of the nonlinear term. The choice of the vector functions and is 
crucial to the success of the method. Mathematically, of course, the sets {\J|j^}'“and 
must be complete in appropriate spaces of functions satisfying (4) plus boundary condi- 
tions and (5) and (6), respectively. From a computation standpoint we would like 
(1) rapid convergence of our expansions of u; (2) miniminn (and ordered) coupling of the 
modes through the time-dependent and viscous operators, for example, banded structure for 
A and B; (3) efficient construction of the matrices A and B; and (4) efficient computa- 
tion of the nonlinear term g. 

In practice, the index n has three components, one for each spatial direction. 

The matrices A and B are diagonal in the homogeneous directions where Fourier series 
may be used but nondiagonal in the direction normal to the wall. For example, in the 
application to pipe flow described below, A and B are banded with the same niimber of 
bands. Thus, implicit treatment of the viscous term is possible at no extra cost. 

Equation (8) is a complete statement of the dynamics. No extra equations are 
required to enforce boundary conditions or continuity, and no fractional time-steps are 
needed. In addition, the number of equivalent grid points in the computation is N/2. 
Thus, only two unknowns per mesh point are required because of the constraint (4) on the 
expansion vectors, allowing considerable savings in memory requirement. 

For pipe flow [16] the boundary condition at the pipe wall r - 1 is u(r ■ 1) ■ 0. 
We assume periodic boundary conditions in the x-direction with period L and write the 
velocity field u as the expansion 

u(r,0,x,t) « 2 ^ 

n,k,i * - 


where each' expansion vector satisfies 

V ■ [^^(r)exp(lkx + iii9)] » 0 


(13) 


and 


X^(l) - 0 . 


(14) 


The weight vectors are given by 

C * Cjjj(r)exp(-ikx - ii6) . (15) 

They are also divergence-free and satisfy the inviscid boundary condition, 

^^(1) . S . 0 . (16) 

. -m r 

[The (k,i) dependence of and is suppressed.] Applying the weighted residual 
method described above we find that for each wave vector (k,Z) 

Aa + ^ Ba - f , (17) 

where 






A “ 

mn 

Jo 

• )S„r dr . 

(18a) 


/•i 



B « 

mn 

1 

• 7 X 7 X Y r dr , 
^n * 

(18b) 
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and 


fm - §0 • \p^ * 2»L ^ sjr dr , 


(18c) 


and where ^ denotes double*Fourler transformation in x and 9. Thus» except for .the 
nonlinear term, the coupling of the equations occurs only through the radial modes as 
discussed above. 


We find that the sequence of expansion vectors *2$o • * •*^n*^n* 
defined in the following satisfies the abovfe requirements. For n > 0 the (r,0,x) com- 
ponents are 



±lkq 


iti 

'n 

Z±i 


kq, 

1 d , 4 1 * 1 / 

*7dF > - 7 ’n 


(k ^ 0) 


(19) 


where 


and (y) Is the shifted Jacobi polynomial [17], 


- rl*’l(l - 


(r^) (n > 0) 


gn^^y) - P^°’'^^2y - 1) , 


( 20 ) 


( 21 ) 


satisfying the orthogonality condition. 


r 


4 (4) / s (4) , „4. 

y 8„ (y)g„ (y)dy - 


( 22 ) 


For n ■ -1 and k ^ 0, we use 


x-i ■ i-i + '^-1 


(23) 


with given by (19) and 


q!i “ rl*’^l - r^) . 


.(24) 


For the case k - 0, the above expansion vectors clearly are not complete and must be 
replaced by an alternative set. A convenient choice Is given by (n > 0) , 


^n 



+ 



(25) 


It Is a simple matter to verify that the expansion vectors defined above yield the 
correct behavior of u(r,k,Z) as r 0; for example, if 1 > 0, 
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- . £-1 
U0 ^ iyv . 




where y 3 ere complex constants. 

The completeness of the expansion Is readily demonstrated. For example, if k 0 
and £ > 0 we find from (19) that 1/2 (u^. t iu0> have the complete representation (and, 
thus, so do and uq, separately), 

■| (u^ t iUg) - ±ika-ir^”^(l - r^) i ik ^n^n"^ * 

nio 

where the first term on the right-hand side of each equation is used to satisfy a nonzero 

value of du0 /dr I that is, (l/ 2 )du 0 /dr | ■ *2ka^j, and the remaining sums are com- 
plete for complex functions having a double zero at the wall and satisfying 

Uy ± iu0 ->■ r^-^ as r 0. The u^ component must satisfy the continuity constraint, 

i£u0 - , 

iku^+— + (ru^) -0 . (27) 


which it does. 


The corresponding weight vectors are essentially the curl of the More spe- 

cifically, if k 0 the weight vectors may be expressed as 



As a result, the (*f) vectors are uncoupled from the (-) vectors. The resulting matrices 
A and B are nonadiagonal, except for an additional nonzero row and a column owing to the 
vector . The limited bandwidth of the viscous matrix B results from the particular 
choice of the polynomials g^^^(y) given by (21). In particular, the Laplacian operator 
in the (r,9) plane is equivalent to a tridiagonal matrix In the following sense: 

7*[q‘(r)exp(lie)] " ^7 ^ ^ q^(r)exp(U0) 


(O, 2 


[b g 
'• n^n- 


1 ( 1 ) 




(r^) + d|^g^^|(r2)]exp(U6) 


In (16] we applied the above technique to the problem of determining the time eigen- 
values for linearized flow in a pipe and demonstrated exponential convergence of the, 
expansion. The results of fully nonlinear calculations of axisymraetric structures are 
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shown In figures 5-7* The evolution of a corotating vortex structure is depicted in 
figure 5, in which azimuthal vorticlty contours are shown* The term "corotating*’ indi- 
cates that the vorticlty in the structure has the same sign as the mean flow vorticlty* 

A counterrotating vortex structure Is shown during various stages of its development in 
figure 6. The counterrotating structure has many features in common with the experi- 
mentally observed turbulent puff [18] — sharp trailing edge, cone-shaped leading edge, 
translation speed -U — but, unlike the experimental puff, which is self-sustaining, it 
decays in time. Clearly, three-dimensional simulations, presently under way, are 
required to capture this intriguing transitional flow structure on the computer. Never- 
theless, the corotating and counterrotating vortex structures depicted appear to be 
important flow structures in more general axisymmetric pipe flow. As shown in figure 7 , 
a random axisymmetric initial field evolves into a flow that is dominated by structures 
of these two types. 

The above method has been applied successfully to channel flow and to flow between 
two concentric cylinders [19]. It seems clear that flow inside a spherical chamber and 
between concentric spheres could be attacked by the same technique, as well as flow in a 
curved pipe. External flows, such as flow past a cylinder or sphere, promise to be more 
challenging because the appropriate behavior at Infinity must be Incorporated into the 
expansion functions and because a rapidly convergent representation of the boundary-layer 
and wake regimes must be provided simultaneously. 

In simple cases, the method can be generalized to the situation where walls are . 
present in more Chan one coordinate direction. For example, in channel flow with 
square (y-z) cross section we might use the expansion vectors. 



'±q;(y)q;(z) - q;(y)q;(z)> 

+ikq^(y)q^(z) 


ikx 




where qj^(s) ■ (1 - s^)^Pj^(s), and is a Legendre polynomial leading to block-banded 

d/dt and viscous matrices with each component matrix banded. For completely arbitrary 
geometries, it appears chat expansion In global basis functions satisfying continuity .and 
the boundary conditions would be cumbersome. In Chat case, using finite-element expan- 
sions with discontintiities in some derivative might be the way to proceed [20]. 


LOW-ORDER EXPANSIONS OF A TURBULENT FIELD 


In the above examples, solutions of the Navler-Stokes equations representing turbu- 
lent flows are computed by essentially tracking the evolution of a point in the (N^) 
dimensional state-space of the coefficients where N is =100 for low Reynolds numbers. 
Can we, by cleverly choosing our basis functions, reproduce the same flow (viz., any 
desired statistic of the flow to within some accuracy) in a much lower dimensional phase 
space? In other words, is the attractor in the 0(10®) space of much lower dimension? 

The answer is probably yes, but we do not yet know how to find these basis functions 
except for some special cases. For example, if the vorticity is concentrated into a 
small region of physical space we achieve a significant reduction in the dimension of , 
coefficient space by tracking vorticity in Lagrangian coordinates [21]. For example, 
one can reproduce the oscillation of a slightly elliptical vortex tube with less than a 

hundred points around the tube. Yet the vortex has a rich energy spectrum, as shown in 

figure 8. Figure 9 shows the three-dimensional wake of a sphere represented by 
^500 points in 10-15 computational vortex filaments [22]. 

If the structure within the vortex tube is an important aspect of the dynamics, then 

several computational filaments must be used to represent the physical vortex tube* 

Ashurst [23] has demonstrated the instability of a vortex ring using this technique. The 
instability of a vortex tube to localized twisting of the vortex lines representing vary- 
ing axial flow within the tube is shown in figure 10. See Nakamura et al. [24] for more 
examples in this area, including examples of vortex breakdown. Chorin [25] has studied 
the statistics of high-Reynolds -number turbulent flow using similar methods. 
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In the general setting, where vorticlty fills a significant portion of the domain, 
Lagrangian vortex methods lose their advantage, requiring as many points as traditional 
Eulerian methods « In this case, the large-eddy simulation (LES) technique [26] is an 
attractive alternative. In LES, one solves the space-filtered Navier-Stokes equations on 
a relatively coarse mesh for the evolution of the large eddies of the turbulent flow. 

The equations contain subgrid-scale stress terms that must be modeled. The method has 
proved quite successful in a variety of applications (see, for example, the LES technique 
applied to channel flow by Moin and Kim [27]) and one can save a large fraction of the 
grid points required to do a direct simulation of the same flow. Even so, it appears 
that one and a half to two decades of scale in each direction (especially if walls are 
present), or a total of half a million points, are needed to represent the important 
energy-containing eddies. 


Clearly it would be advantageous if we could represent, say, turbulent channel flow 
with approximately 100 interacting flow structures, each structure depending parametri- 
cally on a few time-dependent coefficients. The dynamics would be given by coupled non- 
linear ODEs for the coefficients. As pointed out by Saffman [28], the idea of building a 
turbulent field out of a collection of predetermined or organized structures is not new. 
Forty years ago Synge and Lin [29] proposed a representation of isotropic turbulence by a 
(random) superposition of Hill spherical vortices. But this idea was used only as a 
kinematic description of turbulence. For it to become a full dynamical description of an 
evolving turbulent flow, one needs equations of motion for the parameters of the vor- 
tices — position, orientation, size, strength — that accurately represent the dynamics 
of the Navier-Stokes equations. Even if these equations of motion were derived, it 
appears doubtful that only a few hundred such vortices would be required to accurately 
represent a general turbulent flow. 


An analogous situation occurs for the one-dimensional Burgers' equation where exact 
dynamical equations are known for the "pole” decomposition of the solution field [30]. 

In this case, the solution to 


3u 

3t 


TTT + u 


9u 

3x 



(30) 


can be represented by 


u(x,t) - -2v 2 


N 

I 

n*i ‘‘ “n 

where the Zj^ are complex and satisfy the ODEs 


(t) • 


dz„ N . 1 

^ - _2v r; — i — . 

dt w 2_ - Z- 


n m 

mi^n 


(31) 


(32) 


For real solutions u(x,t) the z^ must come in complex conjugate pairs. Thus, each 
flow structure is given by 


-4v(x - Re z) 

(x - Re z)^ + (Im z)^ 


(33) 


which depends on the complex parameter z. However, it appears that the number of struc- 
tures (33) required to reproduce an arbitrary initial field is proportional to the ratio 
of the largest to smallest scale of that field. ^ So although the pole decomposition in 
one-dimensional Burgers' flow has some of the features we are looking for in our flow- 
structure representation of three-dimensional Navier-Stokes flow — exact or approximate 
equivalence to a system of ODEs for the structure parameters — it has the undesirable 
feature that we wish to avoid in three dimensions — the large number of structures 
required to represent an arbitrary flow field. 
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SUMMARY 


The computer simulation of turbulent flows Is a challenging numerical problem. 
Because the three-dimensional turbulent flow field has a wide range of scales with . 
chaotic behavior, one looks for rapidly convergent spatial representations of these 
fields to minimize the number of ODEs in the time reqtilred. In this context, spectral 
expansions appear to be a powerful tool but special attention must be paid to minimize 
the work caused by coupling of the resultant ODEs. A number of low-Reyno Ids -number homo- 
geneous flows have been computed, using meshes with up to 128® points; they have produced 
significant results for understanding turbulence and important information for turbulence 
modelers. We expect results from the direct simulation of wall-bounded turbulent flows 
in the near future, using spectral expansions. In this case, further difficulties arise 
because of the presence of the no-slip condition and the fact that there is a minimum 
Reynolds number, of the order of a few thousand, to sustain turbulence. For these flows, 
expansion of the velocity field in terms of diver gence-free vectors satisfying the , 
no-slip condition appears particularly promising. 

In any event, the direct simulation technique will continue to strain computational 
resources and thereby remain a research tool. More practical low-order expansion tech- 
niques, such as large-eddy simulation and vortex methods, have already made significant 
contributions and will continue to be improved as modeling becomes more sophisticated. 

It is hoped that eventxially modeling will evolve such that only a few hmdred computa- 
tional flow structures are required to represent a complex turbulent flow. 
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TABLE 1.- MESH-POINT REQUIREMENTS 


Reynolds 

ntanber 

Kolmogorov length 

Number of 
mesh points 
(A - 2n) 

N - lOir /D\® 
4 VA/ 

n/D 

n+ 

wall units 

5 X 10* 

0.0045 

1.6 

1.0 X 10^ 

1 X lO** 

0.0028 

1.8 

4.3 X lo’ 

5 X 10" 

0.00093 

2.4 

1.2 X 10® 

1 X 10® 

0.00058 

2.8 

5.1 X 10® 

5 X 10® 

0.00019 

3.8 

1.5 X 10“ 



Figure 1.- One-dlniensional velocity spectra in a pipe. Re « 500,000; r* 
from pipe wall. Dissipation peak is approximate location of max{kil 
Kolmogorov wave number (from Laufer [7]). 


is distance 
i(ki)}; kr^ is 
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E(K» 

EKK) 

E2(K) 

E3(K) ‘ 



CASE 

SHEAR 

RATE VISCOSITY 

o 

BSH9 

28 

0.007 

Q 

BSH10 

28 

0.014 

o 

BSH11 

57 

0.014 

A 

BSH12 

20 

0.005 

♦ 

TC 





- CHAMPAGNE et al. 





Figure 2.- Three-dimeiiiilonal energy spectra for homogeneous shear 
turbulence; 1*^^ ^ is the integral scale in the streamwise direc- 
tion; experiments cited: TC, S. Tavoularis and S. Corrsin» 

J. Fluid Mech. , Vol. 104, pp. 311-347, 1981; F. H. Champagne, 

V. G. Harris, and S. Corrsin, J. Fluid Mech., Vol. 41, 
pp 81-139, 1970 (from Rogallo [121). 


0 2 4 6 8 10 12 14 16 18 20 

TIME. $t 

Figure 3.- The tendency toward equilibrium 
of the turbulence and shear time scales 
for homogeneous shear turbulence (from 
Rogallo [12]). 




(a) Correlation between tensor <(>ij and tensor bij In homogeneous shear 
turbulence (from Rogallo [12]). 



(b) Comparison of measured (4>ij) and modeled (Bbij) terms (from Rogallo [12]). 



(c) Variation of model coefficient with Reynolds number, Rj^ « (from Rogallo [12]). 

Figure 4.- Modeling of the tensor 
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Figure 5,- Evolution of a corotating vortex in pipe flow. 
Contours of azimuthal vorticity; radial scale expanded 
by a factor of 10; Re =* 2200; mean flow is from left 
to right. 



Figure 6.- Evolution of a counterrotating vortex in 
pipe flow. 











Figure 8,^ Three-dimensional energy spectrum of a ring 
vortex; R Is ring radius; F Is circulation j vortex: 
core radius is zero. 


Figure 7*— Evolution of a random axisynm^trlc field 
in pipe flow. * ' 






tu JD » 5.0 7.5 

Figure 9.- Vortex simulation of flow past a sphere with an impulsive start as observed 
from two directions perpendicular to one another (from Leonard [22]). 
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tv 
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Figure 10.- Vortex simulation of the instability of a tube of vorticity because of local- 
ized twisting of the vortex lines. 
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